using CSV,DataFrames,Plots,Statistics,Distributions,Measures,StatsPlots,ColorSchemes
gr();
include("auxilliary.jl"); include("parameters.jl"); include("flow.jl"); include("abcmc.jl");
ENV["COLUMNS"]=200;
# Indices of the ABC files that store the ABC parameter and trajectory files
frt = 1; lst = 10;
# threshold for ABC RMSE error cutoff, given as a top quantile (eg 0.05 is top 5%)
qℓ = 0.1;
# Plot formatting
# Finding that to get right output dimensions and font sizes for gr have to manually
# compute scaling factors based on an observed print output
λplot = 100*7.5/10.42; λfont = 12/18; #scaling factors
# Final output two panel figure dimensions
figwdth = 7.45; fighght = 3.25; # inches
# Final output figure fontsizes
figtickfontsize = 12; figguidefontsize = 14; figtitlefontsize = 16;
lnopt = [:solid, :dash, :dot, :dashdot, :dashdotdot];
mkopt = [:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross,
:xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon,
:heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline,
:+, :x];
df = CSV.read("ABCsmp$frt.csv",DataFrame,header=false);
for id=frt+1:lst
dftemp = CSV.read("ABCsmp$id.csv",DataFrame,header=false);
df = hcat(df,dftemp[:,2:end],makeunique=true);
end
# Dataframe of trajectories
dftrj = CSV.read("ABCtrj$frt.csv",DataFrame,header=false);
for id=frt+1:lst
dftemp = CSV.read("ABCtrj$id.csv",DataFrame,header=false);
dftrj = hcat(dftrj,dftemp[:,2:end],makeunique=true);
end
ncols = ncol(df)-1; nℓcols = qℓ*ncols;
println("The total number of abc sampling before conditioning: $ncols")
println("The total number of abc sampling after conditioning with top $qℓ quantile: $nℓcols")
The total number of abc sampling before conditioning: 5000 The total number of abc sampling after conditioning with top 0.1 quantile: 500.0
prmrg,prmvary=abcdata();prmvary[:ℓerr]=true;prmvary[:βθ]=true;
println("Varied parameters:")
for key in keys(prmvary)
if prmvary[key]
println(key)
end
end
Varied parameters: ρ ℓerr βθ αeff βα γα γθ
# Create dictionary of row position of key parameters;
mykeys = [key for key in keys(prmvary)]; mykeys=vcat(mykeys,[:ℓerr,:βθ])
pnt=Dict{Symbol,Int64}()
for key in mykeys
pos = 1;
while df[pos,1]!=string(key)
pos+=1;
end
pnt[key]=pos;
end;
# Truncate down to trajectories that didn't fail to converge at their numerical tolerances
function Base.isnan(s::Union{Symbol,AbstractString})
return false
end
dfnan = .!(isnan.(df[pnt[:ℓerr]:pnt[:ℓerr],:])); flagnan = [dfnan[1,i] for i=1:ncol(dfnan)];
df = df[:,flagnan]; dftrj = dftrj[:,flagnan];
print("The number of sims that failed to converge is "); println(sum(.!(flagnan)))
The number of sims that failed to converge is 1
# Extract marginal values conditioned on threshold
marginals = Dict{Symbol,Vector{Float64}}();
for key in keys(prmvary)
if prmvary[key]
marginals[key] = [v for v in df[pnt[key],2:end]];
end
end
flag = marginals[:ℓerr] .<= quantile(marginals[:ℓerr],qℓ);
for key in keys(marginals)
marginals[key] = marginals[key][flag]
end;
# Extract trajectories conditioned on threshold
dftrjabc = dftrj[:,2:end];
dftrjabc = dftrjabc[:,flag];
dftrjabc = hcat(DataFrame("day"=>convert(Vector,0:(nrow(dftrjabc)-1))),dftrjabc);
last(dftrjabc,7)
7 rows × 501 columns (omitted printing of 481 columns)
| day | Column3 | Column14 | Column17 | Column22 | Column30 | Column34 | Column39 | Column57 | Column61 | Column63 | Column66 | Column72 | Column85 | Column96 | Column105 | Column106 | Column107 | Column117 | Column130 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 55 | 8287.16 | 6054.01 | 8596.8 | 8241.08 | 8778.35 | 9934.67 | 9871.39 | 8675.37 | 8851.32 | 9589.75 | 9502.8 | 8200.3 | 7159.61 | 9212.67 | 8619.24 | 7965.41 | 9022.96 | 9555.09 | 8301.21 |
| 2 | 56 | 8276.62 | 6001.05 | 8493.72 | 8212.31 | 8766.59 | 10019.4 | 9918.44 | 8669.85 | 8775.63 | 9636.07 | 9557.39 | 8179.52 | 6934.02 | 9143.5 | 8614.97 | 7962.06 | 9077.26 | 9594.97 | 8123.06 |
| 3 | 57 | 8263.3 | 5950.36 | 8382.95 | 8171.35 | 8738.56 | 10095.1 | 9935.27 | 8655.34 | 8703.01 | 9655.54 | 9609.01 | 8149.9 | 6714.52 | 9062.82 | 8606.68 | 7956.75 | 9115.63 | 9612.36 | 7958.81 |
| 4 | 58 | 8252.61 | 5901.46 | 8253.94 | 8123.48 | 8700.8 | 10165.9 | 9917.77 | 8635.0 | 8617.94 | 9647.61 | 9661.05 | 8112.4 | 6486.68 | 8960.48 | 8597.13 | 7951.05 | 9125.09 | 9601.81 | 7800.75 |
| 5 | 59 | 8248.46 | 5853.46 | 8099.61 | 8074.53 | 8660.98 | 10236.0 | 9867.93 | 8612.29 | 8507.5 | 9616.81 | 9715.97 | 8068.97 | 6240.1 | 8829.81 | 8588.49 | 7946.12 | 9098.05 | 9563.43 | 7640.3 |
| 6 | 60 | 8252.28 | 5805.3 | 7917.44 | 8029.44 | 8625.77 | 10308.9 | 9793.02 | 8590.14 | 8363.58 | 9573.7 | 9774.57 | 8021.95 | 5969.88 | 8668.71 | 8581.85 | 7942.52 | 9033.44 | 9494.61 | 7469.87 |
| 7 | 61 | 8263.02 | 5756.22 | 7709.59 | 7991.12 | 8599.23 | 10386.9 | 9703.77 | 8570.38 | 8184.02 | 9533.22 | 9836.17 | 7973.57 | 5677.03 | 8479.69 | 8577.07 | 7940.11 | 8936.6 | 9400.46 | 7284.23 |
histogram(marginals[:ℓerr],title="ℓerr",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:ρ],title="ρ",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:αeff],title="αeff",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:βα],title="βα",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:βθ],title="βθ",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:γα],title="γα",labels="",size=(500,250),normalize=:pdf)
histogram(marginals[:γθ],title="γθ",labels="",size=(500,250),normalize=:pdf)
histogram2d(marginals[:γθ],marginals[:γα],labels="",xlabel="γθ",ylabel="γα",size=(500,250))
βμ = [mean(Weibull(marginals[:βα][i],marginals[:βθ][i])) for i=1:length(marginals[:βα])];
γμ = [mean(Weibull(marginals[:γα][i],marginals[:γθ][i])) for i=1:length(marginals[:γα])];
βσ = [std(Weibull(marginals[:βα][i],marginals[:βθ][i])) for i=1:length(marginals[:βα])];
γσ = [std(Weibull(marginals[:γα][i],marginals[:γθ][i])) for i=1:length(marginals[:γα])];
histogram2d(βμ,γμ,xlabel="Length of infectious period mean",ylabel="Length of recovery period mean",size=(500,325))
plt = palette(:default);
pltβγμ = histogram(βμ,normalize=:pdf,labels="",linecolor=:white,alpha=0.5,linewidth=0);
histogram!(γμ,normalize=:pdf,labels="",linecolor=:white,alpha=0.5,linewidth=0)
density!(βμ,c=plt[1],linewidth=3,linestyle=:solid,labels="contact interval β");
density!(γμ,c=plt[2],linewidth=3,linestyle=:dash,labels="infectious period γ");
plot!(xlabel="mean",ylabel="pdf",xlims=(12.25,15.0),
size=(figwdth/2*λplot,fighght*λplot),
tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
savefig("postmean.pdf")
# Master figure for all posteriors
# Ex of how to use blank plots to center and vary number of subplots in a plot
#pltblank = plot(legend=false,grid=false,foreground_color_subplot=:white);
#lay = @layout [_{0.165w} a{0.335w} b{0.335w} _{0.165w};c d e];
#plot(pltblank,plt1,plt2,pltblank,plt3,plt4,plt5,layout = lay,size=(figwdth*λplot,2*fighght*λplot))
plt1 = histogram(marginals[:ρ],title="ρ",labels="",size=(500,250),normalize=:pdf,ylabel="pdf",xlabel="",
xticks=0.0:0.025:0.05,linecolor=:white,linewidth=0)
plt2 = histogram(βμ,title="βμ",labels="",size=(500,250),normalize=:pdf,xlabel="",
xticks=12.4:0.3:13.0,linecolor=:white,linewidth=0)
plt3 = histogram(βσ,title="βσ",labels="",size=(500,250),normalize=:pdf,xlabel="",linecolor=:white,linewidth=0)
plt4 = histogram(marginals[:αeff],title="αeff",labels="",size=(500,250),normalize=:pdf,ylabel="pdf",xlabel="value",
linecolor=:white,linewidth=0)
plt5 = histogram(γμ,title="γμ",labels="",size=(500,250),normalize=:pdf,xlabel="value",
xticks=12.5:1.0:15.0,linecolor=:white,linewidth=0)
plt6 = histogram(γσ,title="γσ",labels="",size=(500,250),normalize=:pdf,xlabel="value",
xticks=1.8:0.4:3.0,linecolor=:white,linewidth=0)
plot!(plt1,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
plot!(plt2,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
plot!(plt3,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
plot!(plt4,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
plot!(plt5,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
plot!(plt6,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
lay = @layout [a b c;d e f]
plot(plt1,plt2,plt3,plt4,plt5,plt6,layout = lay,size=(figwdth*λplot,2*fighght*λplot))
savefig("postprm.pdf");
The $R_0$ value is given by $R_0 = \int^\infty_0 S_\gamma(t)\beta(t)dt$, where $S_\gamma(t)$ is the survival function of $\gamma$ and $\beta$ is the contact interval hazard. For a Weibull distribution, the survival function is $S_\gamma(t) = e^{-(t/\gamma_\theta)^{\gamma_\alpha}}$. Combining with the definition of $\beta$ we obtain
$R_0 = \int^\infty_0 e^{-(t/\gamma_\theta)^{\gamma_\alpha}}\frac{\beta_\alpha}{\beta_\theta}\left(\frac{t}{\beta_\theta}\right)^{\beta_\alpha-1}dt.$
We now estimate this integral by truncating it to a finite domain with error at most tol, and then use a trapezoidal approximation with npts many sample points.
function R₀(β_α::Real,β_θ::Real,γ_α::Real,γ_θ::Real;tol::Real=1e-6,npts::Integer=20001)
# Find L₀ where tail of integral is at most tol when the tail is bounded by exp(-s/2)
# You want I<tol*β_θ^β_α*γ_α/β_α/γ_θ^β_α
# d/dx F = exp(-s/2) for F = -2 exp(-s/2); so I = 2 exp(-s/2)
L₀ = -2*log(0.5*tol*β_θ^β_α*γ_α/β_α/γ_θ^β_α);
L₀ = L₀>1 ? L₀ : 1.0;
# Find where in [⋅>2r], the exponential is above the poly and may dominate with exp(-s/2)
r = β_α/γ_α-1;
L₀ = L₀ > 2*r ? L₀ : 2*r;
while r*log(L₀) > L₀/2
L₀ *= 2;
end
# Compute definite integral on ∫^L₀_0 by trapezoidal rule
# Note since u^(β_α/γ_α-1) can blow up at 0 in an integrable fashion
# we truncate the origin sample point so an inessential singularity
# does not cause problems
xs = LinRange(L₀/npts,L₀,npts-1); Δ = xs[2]-xs[1];
ys = β_α*γ_θ^β_α/β_θ^β_α/γ_α*exp.(-xs).*(xs.^r);
val = [Δ*(ys[i]+ys[i-1])/2 for i=2:length(xs)];
return sum(val)
end;
marginals[:R₀] = [R₀(marginals[:βα][i],marginals[:βθ][i],
marginals[:γα][i],marginals[:γθ][i]) for i=1:length(marginals[:βα])];
pltR₀ = histogram(marginals[:R₀],normalize=:pdf,labels="",size=(450,250),bins=15,linecolor=:white,linewidth=0,
alpha=0.5)
flag = marginals[:R₀].==Inf; marginals[:R₀]=marginals[:R₀][.!flag];
println("The mean R₀: ",sum(marginals[:R₀])/length(marginals[:R₀]))
println("A 95% CI for R₀: ",[quantile(marginals[:R₀],0.025),quantile(marginals[:R₀],0.975)])
density!(marginals[:R₀],c=plt[1],linewidth=3,labels="")
plot!(xlabel="R₀ value",ylabel="pdf",xticks=0.75:0.5:2.25,
size=(figwdth/2*λplot,fighght*λplot),
tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
The mean R₀: 1.2803727577313408 A 95% CI for R₀: [0.8880207220940958, 1.944244680190886]
savefig("postR0.pdf")
lay = @layout [a b];
plot!(pltR₀,ylabel="");
plot(pltβγμ,pltR₀,layout=lay,
size=(figwdth*λplot,fighght*λplot),
tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
savefig("mst_meanR0.pdf");
dftemp = DataFrame(:βα=>quantile(marginals[:βα],[0.0275,0.975]),
:βθ=>quantile(marginals[:βθ],[0.0275,0.975]),
:γα=>quantile(marginals[:γα],[0.0275,0.975]),
:γθ=>quantile(marginals[:γθ],[0.0275,0.975]),
:ρ=>quantile(marginals[:ρ],[0.0275,0.975]),
:αeff=>quantile(marginals[:αeff],[0.0275,0.975]),
:βμ=>quantile(βμ,[0.0275,0.975]),
:βσ=>quantile(βσ,[0.0275,0.975]),
:γμ=>quantile(γμ,[0.0275,0.975]),
:γσ=>quantile(γσ,[0.0275,0.975])
)
2 rows × 10 columns
| βα | βθ | γα | γθ | ρ | αeff | βμ | βσ | γμ | γσ | |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 5.05714 | 13.5036 | 4.97674 | 13.4295 | 0.00792989 | 0.706579 | 12.4068 | 1.61697 | 12.5304 | 1.93351 |
| 2 | 9.54716 | 13.5076 | 7.9382 | 15.0333 | 0.0495332 | 0.992436 | 12.8249 | 2.81312 | 14.1182 | 2.91921 |
μs = Dict{Symbol,Float64}(); σs = Dict{Symbol,Float64}();
for key in keys(marginals)
μs[key] = sum(marginals[key])/length(marginals[key]);
σs[key] = √( sum((marginals[key].-μs[key]).^2)/length(marginals[key]) );
end
n=length(keys(marginals));
PC = Matrix{Float64}(undef,n,n); id₁=0;
for key₁ in keys(marginals)
id₁+=1;
id₂=0;
for key₂ in keys(marginals)
id₂+=1;
PC[id₁,id₂] = sum( (marginals[key₁].-μs[key₁]).*(marginals[key₂].-μs[key₂]) )/(σs[key₁]*σs[key₂]);
PC[id₁,id₂] *= 1/length(marginals[key₁]);
end
end;
dfpc = DataFrame("prm"=>[key for key in keys(marginals)]);
pos = 0;
for key in keys(marginals)
pos += 1;
dfpc[!,string(key)] = convert(Vector,PC[:,pos]);
end
println("Pearson correlation coefficients:")
dfpc
Pearson correlation coefficients:
8 rows × 9 columns
| prm | R₀ | ρ | ℓerr | βθ | αeff | βα | γα | γθ | |
|---|---|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | R₀ | 1.0 | 0.593828 | 0.709237 | 0.453027 | 0.0258399 | 0.453027 | -0.106132 | 0.792204 |
| 2 | ρ | 0.593828 | 1.0 | 0.0189204 | 0.349604 | 0.0488572 | 0.349573 | -0.069226 | 0.44916 |
| 3 | ℓerr | 0.709237 | 0.0189204 | 1.0 | 0.383182 | -0.00831572 | 0.383154 | -0.0803331 | 0.592233 |
| 4 | βθ | 0.453027 | 0.349604 | 0.383182 | 1.0 | 0.0756349 | 1.0 | 0.124744 | -0.067327 |
| 5 | αeff | 0.0258399 | 0.0488572 | -0.00831572 | 0.0756349 | 1.0 | 0.0756013 | 0.0805348 | -0.00315127 |
| 6 | βα | 0.453027 | 0.349573 | 0.383154 | 1.0 | 0.0756013 | 1.0 | 0.12476 | -0.0673385 |
| 7 | γα | -0.106132 | -0.069226 | -0.0803331 | 0.124744 | 0.0805348 | 0.12476 | 1.0 | 0.0897513 |
| 8 | γθ | 0.792204 | 0.44916 | 0.592233 | -0.067327 | -0.00315127 | -0.0673385 | 0.0897513 | 1.0 |
df[!,1] = Symbol.(df[!,1])
# Find best fit
idℓerr = 1;
while df[idℓerr,1]!=:ℓerr
idℓerr+=1;
end
ℓs = [df[idℓerr,k] for k=2:ncol(df)]
ℓbf = minimum(ℓs);
pos = findfirst(ℓs.==ℓbf) + 1;
# Prepare the data vector
vkeys = df[:,1]; prm,_=rdprm(df[:,pos],vkeys);
data!(prm);
println("The best fit parameters:")
for key in keys(prm)
println(key,": ",prm[key])
end
The best fit parameters: αeff: [0.9698334727742906] nndsmp: [2500.0] αL: [14.0] nnd: [2500.0] dwnsmp: [1.0] βθ: [13.503792220531302] atol: [1.0e-5] Δtmin: [1.0e-5] ℓerr: [746.6170971719679] γα: [6.052065267527739] γθ: [13.617999656279881] fˢη: [3.301535733241945e-5] yˢrg₀: [0.0, 36500.0] yᵛrg₀: [0.0, 0.0] λ: [1.5] βα: [5.238644772185896] yⁱrg: [0.0, 22.0] ρ: [0.04476143286480499] nΔtfail: [NaN] rtol: [0.001] yᵛrg: [0.0, 22.0] yˢrg: [0.0, 36565.0] fⁱη: [0.07407407402502883] yⁱrg₀: [0.0, 14.0] Trg: [0.0, 65.0]
# Run the simulation
ysol,yʳsol = pdesolve(;prm=prm);
Simulation progress: 0.015625/1 Simulation progress: 0.03125/1 Simulation progress: 0.046875/1 Simulation progress: 0.0625/1 Simulation progress: 0.078125/1 Simulation progress: 0.09375/1 Simulation progress: 0.109375/1 Simulation progress: 0.125/1 Simulation progress: 0.140625/1 Simulation progress: 0.15625/1 Simulation progress: 0.171875/1 Simulation progress: 0.1875/1 Simulation progress: 0.203125/1 Simulation progress: 0.21875/1 Simulation progress: 0.234375/1 Simulation progress: 0.25/1 Simulation progress: 0.265625/1 Simulation progress: 0.28125/1 Simulation progress: 0.296875/1 Simulation progress: 0.3125/1 Simulation progress: 0.328125/1 Simulation progress: 0.34375/1 Simulation progress: 0.359375/1 Simulation progress: 0.375/1 Simulation progress: 0.390625/1 Simulation progress: 0.40625/1 Simulation progress: 0.421875/1 Simulation progress: 0.4375/1 Simulation progress: 0.453125/1 Simulation progress: 0.46875/1 Simulation progress: 0.484375/1 Simulation progress: 0.5/1 Simulation progress: 0.515625/1 Simulation progress: 0.53125/1 Simulation progress: 0.546875/1 Simulation progress: 0.5625/1 Simulation progress: 0.578125/1 Simulation progress: 0.59375/1 Simulation progress: 0.609375/1 Simulation progress: 0.625/1 Simulation progress: 0.640625/1 Simulation progress: 0.65625/1 Simulation progress: 0.671875/1 Simulation progress: 0.6875/1 Simulation progress: 0.703125/1 Simulation progress: 0.71875/1 Simulation progress: 0.734375/1 Simulation progress: 0.75/1 Simulation progress: 0.765625/1 Simulation progress: 0.78125/1 Simulation progress: 0.796875/1 Simulation progress: 0.8125/1 Simulation progress: 0.828125/1 Simulation progress: 0.84375/1 Simulation progress: 0.859375/1 Simulation progress: 0.875/1 Simulation progress: 0.890625/1 Simulation progress: 0.90625/1 Simulation progress: 0.921875/1 Simulation progress: 0.9375/1 Simulation progress: 0.953125/1 Simulation progress: 0.96875/1 Simulation progress: 0.984375/1 Simulation progress: 1.0/1 Number of times Δt was too small for refinement: 0
plot(:α;prm=prm)
plot(:β;prm=prm)
plot(:γ;prm=prm)
plot(:Weibull;prm=prm)
plot(:fˢ;prm=prm)
plot(:fⁱ;prm=prm)
plot(ysol)
plot(ysol,yʳsol;prm=prm)
plotbd(ysol;prm=prm)
df_yⁱ = CSV.read("ODH_snipdaily.csv",DataFrame);
first(df_yⁱ,3)
3 rows × 12 columns
| time | 0-9 | 10-19 | 20-29 | 30-39 | 40-49 | 50-59 | 60-69 | 70-79 | 80+ | daily_confirm | cum_confirm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Date | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 2020-11-15 | 559.202 | 578.048 | 1628.5 | 1341.0 | 1335.75 | 1421.5 | 1134.5 | 679.25 | 474.75 | 9152.5 | 3.60476e5 |
| 2 | 2020-11-16 | 561.538 | 580.462 | 1642.2 | 1354.4 | 1344.0 | 1432.4 | 1138.4 | 676.2 | 486.2 | 9215.8 | 365386.0 |
| 3 | 2020-11-17 | 556.457 | 575.21 | 1641.17 | 1368.33 | 1363.33 | 1440.5 | 1154.5 | 685.167 | 499.667 | 9284.33 | 370264.0 |
# Extract error for this trajectory
rms,yⁱ_daily = ℓerr(ysol;prm=prm);
# Compute 95% quantiles
qlow = [quantile(dftrjabc[i,2:end],0.0001) for i=1:nrow(dftrj)];
qhgh = [quantile(dftrjabc[i,2:end],0.9999) for i=1:nrow(dftrj)];
npts = nrow(df_yⁱ);
#plot(df_yⁱ[!,:time],df_yⁱ[!,:daily_confirm],linewidth=1,labels="ODH raw",linestyle=:dash,
# title="ODH Daily Cases",xlabel="date",ylabel="count",size=(400,225));
odhavg = Vector{Float64}(undef,npts);
@inbounds for i=1:npts
Σ=0.0; ct = 0;
for k=-3:3
try
Σ += df_yⁱ[i+k,:daily_confirm];
ct += 1;
catch
true;
end
end
odhavg[i] = Σ/ct;
end
plt = palette(:default);
plot(df_yⁱ[!,:time],yⁱ_daily[1:npts],linewidth=0,labels="",c=plt[1],
ribbon=(yⁱ_daily[1:npts]-qlow,qhgh-yⁱ_daily[1:npts]),
fillalpha=0.25)
plot!(df_yⁱ[!,:time],yⁱ_daily[1:npts],linewidth=3,labels="model",c=plt[1],
#ribbon=(yⁱ_daily[1:npts]-qlow,qhgh-yⁱ_daily[1:npts]),
xlabel="date",ylabel="daily cases")
plot!(df_yⁱ[!,:time],odhavg,linewidth=3,labels="ODH 7 day mvavg",c=plt[2],
linestyle=:dash);
plot!(legend=:bottomright)
plot!(size=(figwdth*0.75*λplot,fighght*λplot),
tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
titlefontsize=round(figtitlefontsize*λfont))
savefig("abcvar.pdf");